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Abstract 

We present a new algorithm that produces a well-spaced superset of points conforming 
to a given input set in any dimension with guaranteed optimal output size. We also provide 
an approximate Delaunay graph on the output points. Our algorithm runs in expected time 
(3(2°( d ) (n log n + m)), where n is the input size, m is the output point set size, and d is the 
ambient dimension. The constants only depend on the desired element quality bounds. 

To gain this new efficiency, the algorithm approximately maintains the Voronoi diagram 
of the current set of points by storing a superset of the Delaunay neighbors of each point. 
£S) By retaining quality of the Voronoi diagram and avoiding the storage of the full Voronoi 

diagram, a simple exponential dependence on d is obtained in the running time. Thus, if 
one only wants the approximate neighbors structure of a refined Delaunay mesh conforming 
to a set of input points, the algorithm will return a size 2°^m graph in 2°^(n\ogn + m) 
expected time. If m is superlinear in n, then we can produce a hierarchically well-spaced 
c/3 superset of size 2°^n in 2°^n\ogn expected time. 

^h 1 Introduction 

> 

Delaunay meshes are used extensively in graphics and scientific computing, and more recently, 
ly-j they have been applied in higher dimensions for data analysis [31 [10]. However, there are 

serious difficulties in constructing meshes in more than three dimensions. Perhaps the two most 
pressing problems are the combinatorial blowup in size complexity and the numerical instability 
of the predicates used in the construction. We present a new algorithm that alleviates both of 
these problems. It works by only storing an approximation to the Delaunay graph. This has 
the simultaneous effect of eliminating much of the complexity of the mesh itself while at the 
same time eliminating the costly and unstable predicate computations. It also introduces a new 
walk-based point location method that can locate the nearest neighbor of a point to be inserted 
in constant time after 2 ( 'niogn time for preprocessing. 

The specific class of meshes we consider arises from Delaunay or Voronoi refinement. The 
input is a set of n points P C M and the output is a well- spaced superset M of size m. 
Delaunay refinement meshes have an extensive theory providing guarantees including optimal 
runtime, optimal output size, conformity to input boundaries, and even the removal of a class 
of poorly-shaped simplices called slivers. These methods are widely used in practice and form 
the basis for several popular meshing codes such as Triangle [18J, TetGen |19j . and the CGAL 
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3D Mesh Generation package [lj . Still, these codes have not been extended to higher than three 
dimensions. 

Delaunay refinement is inherently incremental. At every point in time, the data structure 
represents a valid Delaunay triangulation. This makes several things easier, but the convenience 
comes at a cost. The numerical computation of linear predicates often gives unstable solutions 
requiring exact arithmetic or other machinations to resolve. A significant amount of research 
has gone into more robust and efficient predicates and yet, for Delaunay refinement, only a 
small fraction of these computations will determine the structure of the output triangulation. 
However, the intermediate state must be exact to guarantee that the output is correct. Similarly, 
representing all of the simplices is a waste as only a fraction of these simplices appear in the 
output. 

We propose to perform a kind of Voronoi refinement that stores only a sparse graph as its 
intermediate state. This graph is an approximation to the edges of the Delaunay triangula- 
tion. Two new challenges arise when working with the approximate Delaunay graph. The first 
challenge is that the circumcenters of Delaunay simplices are usually used to refine the mesh. 
Without simplices, one must look elsewhere for good places to insert Steiner points. The second 
challenge is that the intermediate meshes are often used as point location data structures for 
the input points that have not been inserted yet. Without the clear decomposition of space 
given by the Delaunay triangulation or the Voronoi diagram, one must find some other way to 
organize the uninserted points geometrically. 

For Steiner points, we show how approximate linear programming can give approximate 
Delaunay circumcenters without ever finding a single Delaunay simplex. For point location, we 
show how preprocessing the points in a special ordering allows for a constant time per point 
scheme. Thus, after 2 ( 'nlogn preprocessing time, the algorithm will run in output-sensitive 
linear time. 

Our Contribution We present a new algorithm that computes an asymptotically optimal 
well-spaced superset of a set of points P. It has all the quality guarantees achieved by Delaunay 
or Voronoi refinement, but never actually constructs the full Delaunay triangulation or Voronoi 
diagram. The expected running time is 2 ' '(nlogn + to), and the output graph has size 
2 ( 'm. If a hierarchically well-spaced point set is sufficient, then the expected running time is 
2°Wnlogn, and the output size 2°^ d >n. 

Related 'work Boissonat et al. [2J proposed to store only the edges of the Delaunay trian- 
gulation as a way to make the construction more tractable in higher dimensions. Still, their 
approach requires enumerating the higher dimensional simplices locally when making updates 
and thus still requires some high dimensional predicates. In order to avoid this work, it seems 
necessary to relax the condition that the exact Delaunay graph be stored. 

Well-spaced points are a standard method in mesh generation. The use of Delaunay re- 
finement to produce well-spaced points dates back to Chew [1]. Ruppert [15] developed the 
optimality theory for such meshes and the generalization of his quality criteria was formalized 
by Miller et al. [13]. The recent book by Cheng et al. [3] gives a thorough treatment of Delaunay 
meshing in M 3 . 

The Sparse Voronoi Refinement algorithm of Hudson et al. [9] changed the perspective 
from Delaunay refinement to Voronoi refinement and was able to prove the first near-optimal 
output-sensitive algorithm with a running time of Od(n log A + to). Here, Od hides constant 
factors that depend only on d. The term depending on the spread A was removed by Miller et 
al. [12] . who devised the NetMesh algorithm, an optimal comparison-based algorithm running 
in Od{n log n + m) time. For d-dimensional quality meshes that include all simplices, these 
constants are necessarily superexponential in d. We avoid these superexponential constants by 
only storing a graph. 



2 Background 

Geometry Basics We assume throughout that all points lie in ti-dimensional Euclidean space. 
We identify the points of this space with vectors in M rf and use d(u, v) to denote the Euclidean 
distance between u and v. Similarly, we write d(u,S) for a point u and a compact set S to 
denote the minimum distance from u to a point of S. The nearest neighbor of u in S is 

NN(u, S) : = argmind(u, v). 
v£S\{u} 

If the nearest neighbor is not unique, we choose NN(u, S) arbitrarily among the possibilities. 
When the underlying set S is clear from context, we will often simply write NN(u). 

Voronoi Diagrams and Delaunay Triangulations The Voronoi diagram, Vorp, of a 
set of points P in M rf is a decomposition of M. d into polyhedral cells called Voronoi cells, where 
each cell is the set of reverse nearest neighbors of a point of P. That is, Vorp(p) = {x £ 
M rf : d(x,p) = d(x,P)}. These cells are closed and form a cell complex, i.e. the intersection of 
two cells is another polytope of one less dimension. The cells of all dimensions of the Voronoi 
diagram are called the faces. The dual of the Voronoi diagram is the Delaunay triangulation, 
Delp. It has a /c-simplex for every (d — /c)-face of the Voronoi diagram. Moreover, the center of 
the ball circumscribing any ci-simplex of Delp is a vertex of Vorp. 

The total number of faces of the Voronoi diagram, or dually, the number of simplices of the 
Delaunay triangulation can be as large as 0(n' ' 2 '). However, if the points are more evenly 
spaced, the complexity drops to 0(n). One view of Voronoi- or Delaunay-based mesh generation 
algorithms is that they seek to find such an arrangement of points for which this complexity 
is as small as possible. However, even in such instances, the dependence on the dimension is 
quite bad. In particular, the complexity is dominated by the higher dimensional simplices. This 
motivates our approach of storing only the edges of Delp. 

Mesh Generation and Well-Spaced Points The inradius of a Voronoi cell Vorp(tt) is the 
radius of the inball, the largest ball centered at u that is contained in Vorp(u). It is denoted 
r(u) and is equal to half the distance from u to its nearest neighbor in P. The outradius of 
Vorp (it), denoted R(u), is the radius of the smallest ball that contains all of the vertices of 
Vorp(ii). The aspect ratio of Vorp(w) is defined as aspectp(tt) = R(u)/r(u). We say that M is 
r-well-spaced if for all u E M, aspect M (u) < r. In such cases, the constant r is sometimes 
called the quality of the mesh Delju- 

For some point sets P, any well-spaced superset M of P will be prohibitively large. To 
combat this problem, Miller et al. defined hierarchically well-spaced points [12]. A set M is 
hierarchically r-well-spaced if there is a rooted tree T such that following hold. 

1. The nodes V of T are subsets of M. 

2. Each set in V is r-well-spaced. 

3. For any two A, B G V, A n B is at most a single point. If A is the parent of B then A n B 
is a single point. If A n B is a point then that point is also in every set in the path from 
A to B in T. 

4. If A is the parent of B in T and p = A n B then B C ball(p, ed(v, A \ {p})) for some small 
constant e that will depend on the details of the algorithm. 

Hierarchically well-spaced points provide a Delaunay refinement analogue of compressed quadtrees. 
Most importantly, for any P C M. d , there exists a hierarchically well-spaced superset of size 
2°W n [J2]. For a wide class of reasonable inputs, the hierarchy is not necessary, in which case 
the output hierarchy will have a single set of linear size. However, for some inputs, no linear-size 
(non-hierarchical) well-spaced superset is possible [17j . 



The Feature Size For any set S, fg maps a point x G M rf to the distance to the second 
nearest point to x in S. For example, if x G S then fs(^) = d(x, NN(x)) since re itself is the 
nearest point in S. This 1-Lipschitz function is known as the feature size with respect to S. 
Many standard guarantees of meshing algorithms are derived from looking at the feature size 
function. For example, to show that an algorithm terminates with an optimal size output, it 
suffices to show that fp(v) < Kfnj{v) for each v G M, where P is the input and M is the 
output. 

Meshing Guarantees There are several desirable guarantees for an algorithm that produces 
a well-spaced superset M of an input set P: 

• Quality Output: aspect m(v) < r for all v G M. 

• Local Sizing: ip(v) < K?m(v) for all v G M and a constant K that depends only on the 
quality r. 

• Size Optimality: \M\ = 2 ^ >\M'\ for any other r-well-spaced superset M' of P. 

• Running Time: Od(n log n + \M\). 

In this paper, we give an algorithm that achieves all of these guarantees. Moreover, the constants 
in the running time are at most 2°^ d > . 

Delaunay and Voronoi Refinement The Delaunay refinement paradigm for mesh gener- 
ation is an incremental algorithm that improves mesh quality by iteratively adding a Steiner 
point at the circumcenter of a simplex whose circumradius to shortest edge ratio is greater than 
some constant. The dual approach, known as Voronoi refinement, instead adds the vertex u 
of Vor(u) farthest from v whenever Vor(u) has aspect greater than some constant r. To avoid 
confusion with the input vertices, we call such a point u the farthest corner of the Voronoi cell 
Vor(w). Recall that vertices of the Voronoi diagram are circumcenters of Delaunay simplices, 
but the quality conditions are slightly different. Although asymptotically the same, Voronoi 
refinement will produce stronger guarantees on the degree of the Delaunay graph and the edge 
length discrepancy of edges incident to a vertex at the cost of doing more refinement. 

The key idea in both methods is to perform only local operations. Then, the structure of 
the current Delaunay triangulation or Voronoi diagram is used to find large empty regions to 
refine. 

Samples, Nets, and Cages A subset S C P is an (e, 5)-sample of P if all points of P are 
within distance e of a point of S and no pair of points of S is separated by a distance less than 5. 
Such samples are also known as Delone sets [5] . The special case of e = 5 is sometimes referred 
to as a (metric) e-net. Note that such samples are not the same as well-spaced points as they 
are restricted to a subset and also cannot adapt their density locally as well-spaced points do. 
In several places in the algorithm we need a net on the (d — l)-sphere Sd—i- We will refer 
to these special nets as cages to avoid confusion with range space nets. We state the following 
folklore result (see Matousek [11] for a proof). 

Lemma 1. For any rj G (0,1] there exists a subset of points C C Sd-i called an n-cage such 
that for all x G S d ^, d(x, C)<r, and \C\ < 0({4/i]) d ). 

We will assume that we have precomputed a fixed 77-cage for sufficiently small rj and refer 
to it as 77— CAGE. Since only one such cage is needed, we assume it is hardcoded into the 
algorithm. 



Permutations and Greed A permutation on a set of points induces an order relation. For a 
set S with an order relation -< and a point v E S, the predecessor radius of v with respect to 
S is defined as X(v, S) := ram u ^ v d(u, v). It is often convenient to express an ordering on a set 
by labeling its points. For example, we write a permutation onPaspi,... ,p n . The prefixes 
of this ordering are the sets P{ := {p±, . . . ,Pi}. A predecessor pairing <j) maps each point pi 
to <j)(pi) = NN(pi, Pt-i) for i > 2. 

The greedy permutation of a set of points P is an ordering that starts with any point p\ 
and then satisfies pi = argmax E pd(p, -Pj-i). This is also sometimes known as a farthest point 
sampling strategy. 

A good approximation to the greedy permutation as well as a predecessor pairing can be 
computed in 2 "nlogn time using an algorithm of Har-Peled and Mendel [8J. It is approximate 
in the sense that each point is approximately the farthest next point up to a factor of 1 — n~ a for 
some a > 0. To simplify the exposition, we assume that we compute a true greedy permutation. 
An approximation will only affect some constant factors. 

Approximate Linear Programming Our algorithm will employ linear programming (LP) 
using, for instance, the ellipsoid method to find large empty regions in the point set corre- 
sponding roughly to the farthest corners of poor quality cells. The running time of the ellipsoid 
method for d dimensions and I constraints is 0(£poly(d) log i), where e is the error bound [7J. 
For the linear programs we care about, it suffices to compute the answer to a fixed constant 
error. Moreover, the number of constraints will always be 2°^ d > . Thus, the time required to 
solve such an LP is also 2 ( >. 

3 The Algorithm 

3.1 Overview of the Algorithm 

The basic structure of the algorithm is as follows: First, one computes a greedy permutation of 
the input points using the algorithm of Har-Peled and Mendel [8] ■ The points are then inserted 
one at a time, yielding an incremental construction of a hierarchically well-spaced set of points. 
Each time an input point is inserted near an existing Steiner vertex q, the vertex q is removed 
in a process we call snapping. After each input point insertion, Steiner vertices are added as 
part of a refinement procedure to retain quality. At each step in the algorithm, we maintain 
an approximate Delaunay graph as well as the inradius and approximate farthest corner of the 
inserted vertices. Like walk-based algorithms for Delaunay triangulation, there is no auxiliary 
data structure required for point location [6j U3], but using greedy permutations allows us to 
prove stronger guarantees. Pseudocode for these routines may be found in Appendix \A\ 

3.2 Data Structures 

The algorithm must maintain a hierarchy of well-spaced points and a sparse graph on each point 
set. The hierarchy is represented as a tree T. The nodes of T are identified with the point sets. 
Recall that in the definition of a hierarchically well-spaced point set, the individual sets may 
intersect at exactly one point, and when this happens, then one set is the ancestor of the other. 

The main data structure used in the algorithm is a sparse graph representing the approx- 
imate Delaunay graph of the point sets in the hierarchy. The neighbors of a point p are 
denoted N{p). The edges of the approximate Delaunay graph contain the true Delaunay edges 
of the points. The extra edges incident to a point p are all nearly Delaunay in that they are 
approximately the same length as the true Delaunay edges incident to p among the points in 
the same set of the hierarchy. 

Rather than a single approximate Delaunay graph, we keep a hierarchy of them, but we 
treat it as a single disconnected graph. However, each point gives a unique handle into one of 



the components; it is the lowest node in the tree containing that point. Hence, there is no need 
to navigate this hierarchy during the algorithm as it is merely an organizational tool. 

Unlike the previous work on hierarchically well-spaced points, the hierarchy grows mono- 
tonically, so, there is no need to maintain the hierarchy under dynamic deletions [12) . As a 
post-process, the collection of approximate Delaunay graphs is combined to produce a single 
approximate Delaunay graph. 

Attached to each vertex of each set in the hierarchy is an inradius r(p), an approximate 
outradius R(p), and an approximate farthest corner. 

3.3 Pre-processing 

The first step of the algorithm is to pre-process the input points by computing an approximate 
greedy permutation p\ , . . . , p n of P. This will be the order in which we insert the input points. 
Let 4> be the corresponding predecessor pairing. 

We start by placing the input points in a bounding domain. This will ensure that we do not 
insert too many points, because otherwise, local refinement rules can generate Steiner points 
to fill all of M. d . We will use a copy of 77— CAGE as the bounding domain scaled to have a 
radius that is a constant times larger than the diameter of P. It is known that for a reasonable 
constants, this will contain the refinement (see [HI Ch. 3] for a detailed discussion). 

3.4 Point Location 

Point location is needed only to identify the nearest neighbor of a point we are trying to insert 
among those vertices that we have already inserted. For a point p, a greedy walk routine 
locates this nearest neighbor of p. We shoot a ray from the predecessor (f){p) to p and keep track 
of the Voronoi cells that the ray passes through. Since we insert the points in the order given 
by the greedy permutation, <p(p) will have been inserted previously. Even though we don't store 
the full Voronoi diagram, the approximate Delaunay graph suffices to find the next Voronoi cell 
that the ray passes through and, thus, the next step in the walk. 

Let v be the vertex of the current position of the walk. The walk begins at v = <j)(p). For 
each step, compute the intersection of the line through p and (j){p) with the bisecting hyperplane 
between q and v for each neighbor q of v. Computing these intersection points only requires 
scalar products rather than the determinant predicates. Among those intersection points for 
which the ray is leaving Vor(t;), choose the intersection that is closest to v. This is the true 
intersection of the boundary of Vor(u) with the ray and indicates the next cell in the walk. The 
walk terminates when the intersection is not in p(p(p), and the last vertex is NN(p). As we show 



in Lemma 14, each greedy walk only takes a constant number of steps. 

The intersection of the halfspaces bounded by orthogonal bisectors of the Delaunay edges 
incident to a vertex is precisely the Voronoi cell of the vertex. Consequently, the correctness of 
this algorithm only requires that the approximate Delaunay graph contains the true Delaunay 
graph as a subgraph. The extra edges cannot change the walk because their orthogonal bisector 
hyperplanes all lie strictly outside the Voronoi cell. 

3.5 Incremental Updates 

The insertion routine adds a new point p to the current data structure. It first finds the nearest 
neighbor of p among those points inserted thus far using the greedy walk procedure described 
above. We then find the neighbors of p in the approximate Delaunay graph by searching locally 
in the graph starting from NN(p). For input points, insertion may cause a nearby Steiner point 
to be deleted or it may cause the creation of a new layer in the hierarchy. These steps are 
described in more detail below. 



Inserting Input Points Let p be an input point to be inserted. Once we have performed 
point location and determined q = NN(p), we perform one of three actions, depending on q and 
the local sizing constant K. 

• If q is another input point and d(p, q) > j^r(q) : we insert p. 

• If g is another input point and d(p,q) < -^r(q), we add a new layer to the hierarchy (see 
below) and add p to the new layer. 

• If q is a Steiner point that has been added during a refinement stage, we snap q to p. 
That is, we delete q and insert p. 

In all cases, we finish by updating the approximate Delaunay graph as described below. 

Adding a New Set to the Hierarchy A new set of the hierarchy is always added as a new 
leaf to the hierarchy tree. This happens only when attempting to add an input point p that is 
too close to another input point q compared to the feature size at q. In such cases, we replace 
the point q with a copy of itself and move q to the newly created set. We then add a new copy 
of rj— CAGE centered at q with radius 2 d(p, q) and add p and q to the new level. This is a new 
connected component of the approximate Delaunay graph and starts as a complete graph. This 
roughly corresponds to restarting the algorithm with just the points p and q, except there is no 
need to recompute the greedy permutation. 

If an input point has nearest predecessor q, then this refers to the copy of q in the new 
set in the hierarchy. This convention makes it possible to ignore the hierarchy when doing 
point location as the greedy permutation guarantees that any changes near q will happen at the 
smaller scale and thus in the newly created layer. 

Finding Approximate Delaunay Neighbors The insertion of p, along with the possible 
removal of q = NN(p) in the case of snapping, alters the current approximate Delaunay graph G. 
To find the neighbors of p in the updated graph, we do a breadth-first search from q in G. The 
search is pruned at any vertex v that cannot be a neighbor of p because either d(p, v) > 2R(v) 
or d(p, v) > 2r"r(p), where r" is an upper bound on the aspect ratio of any vertex at any time 



in the algorithm. The number of neighbors found by this search is at most 2 ™ > (Lemma 10 ) 
They are used to estimate the outradius R(p) after inserting p (as described below). The list of 
new neighbors is then limited to those vertices v found in the search such that d(p, v) < 2R(p), 
as v cannot otherwise be a Delaunay neighbor of p. Only the neighbors of p need to be updated 
by recomputing inradii, outradii, and approximate farthest corners. 

When deleting a point p, we update the approximate Delaunay graph by adding edges 
between all pairs (u, v) of neighbors of p. Then, we recompute the approximate aspect ratios of 
the points and eliminate any edges for which d(u,v) > 2mm{R(u), R(v)}. 

3.6 Approximate Farthest Corner 

For each point p whose approximate Delaunay neighborhood is updated, we update the approx- 
imate farthest corner of Vor(p) as follows. 

Given a unit vector c G S 1 , we can find the corner of the Voronoi cell of p that is farthest 
in the direction of c by solving a linear program. The Voronoi cell of p is precisely the set of 
points satisfying a linear constraint for each q £ M(p) given by the hyperplane bisecting pq. 
Thus, the LP with these constraints and objective function given by the inner product with c 
has a solution on the boundary of Vor(p) that is extremal in the direction of c. 

If c is chosen to be sufficiently close to the unit vector in the direction of the farthest corner 
v , then v maximizes the objective function of our linear program. However, the problem is 
that we do not have any a priori knowledge of the optimal vector c. Thus, we will solve the 
aforementioned problem for several choices of c. In particular, we iterate c over all points in 
rj— CAGE in order to guarantee that some c is close to the direction of v. 



To solve the LP, we can either use an interior point or ellipsoid method, as both would run 
in polynomial time in d and the number of constraints. In our case, the dimension is quite small 
while the number of constraints is exponential in the dimension. Thus, asymptotically we get 
better bounds using the ellipsoid method. Moreover, the aspect ratio bounds make it easy to 
find a sufficiently good starting ellipse. We discuss a scheme based on the ellipsoid method for 
our problem. 

In the ellipsoid method, we need an ellipse that contains the Voronoi region, but in our 
application the aspect ratio of the Voronoi cells is at most a constant r". Thus, for a given point 
p, ball(p, r{p)r") contains the feasible region. Moreover, we shall only require an approximate 
solution, namely, a point x for which c T x is within (1 — e) of the optimum. Thus, for some 
value of t, we add the additional constraint ex > t to our LP. We then search for a feasible 
solution. The trick is to perform a binary search on t £ [r(p), r(p)r"] and solve the LP feasibility 
for different values of t. Since we desire only an approximate solution within a (1 — e) factor, 
we obtain a lower bound on the volume of the feasible region and can bound the number of 
iterations of the ellipsoid method required. 

Each LP feasibility computation will require time 0{m + poly(d)) log(r"/e), where m is the 
number of constraints (in our case, 2°^ d >). The binary search will add an additional multiplica- 
tive factor of log(r"/e). Moreover, repeating the procedure for all c G r/— CAGE will incur an 
additional 2°^ factor. 

3.7 Refinement 

The insertion of a new input point may cause some Voronoi cells to have aspect ratio more than 
r. The refinement phase adds Steiner points in order to recover quality. 

Having computed the approximate farthest corner for all updated vertices, we have an 
approximate aspect ratio for each Voronoi cell. We wish to keep a queue of all Voronoi cells 
whose aspect ratio is greater than r. Since we only know approximate aspect ratios, we instead 
add to the queue all Voronoi cells for which our computed approximate aspect ratio is greater 
than t(1 — e) ( 1 — \\- This is because any Voronoi cell with a true aspect ratio of more than 
r must necessarily satisfy this criterion, as is shown in Corollary [T] 

The refinement stage consists of going through the queue of bad aspect ratio cells and 
inserting the approximate farthest corner of each as a Steiner vertex. If a vertex is added to 
the queue but some other refinement causes it to no longer have bad aspect ratio by the time 
we process it, then we simply do nothing and move on to the next cell in the queue. Each 
insertion requires that we update the approximate Delaunay graph exactly as we did for input 
points. The point location cost for Steiner vertices is constant time per vertex because the 
nearest neighbor of the Steiner point is just the vertex whose cell had bad aspect ratio in the 
first place. 

Refinement can cause new cells to have bad aspect ratio. These are also added to the queue 
after each insertion. Refinement continues until the queue of Voronoi cells is empty. 

3.8 Post-processing 

The hierarchy can be reconciled into a connected approximate Delaunay graph as a post-process 
in one of two ways. A simple, output-sensitive linear-time algorithm is known for extending a 
hierarchically well-spaced point set into a single well-spaced set [12]. The same algorithm can 
be adapted for our algorithm as it only requires the ability to refine Voronoi cells. However, 
since this operation could potentially add a superlinear number of Steiner points, we instead 
opt to complete the approximate Delaunay graph by adding edges only. 

Let A, B be two nodes of the hierarchy tree such that A is the parent of B. Let v be the 
vertex shared by both A and B and let Ga and Gb be the approximate Delaunay graphs of 
A and B respectively. Let C be the vertices in the boundary cage of level B and let D be the 



vertices adjacent to v in Gb- We replace the nodes A and B with a new node whose vertex set 
is A U B and its edge set contains the union of the edges in Ga and Gb as well as the complete 
bipartite graph on (C,D). The new tree is the formed from the old one by contracting the 
edge from A to B. Continuing this way allows us to flatten the entire hierarchy into a single, 
connected approximate Delaunay graph. 

4 Algorithmic Invariants 

The analysis follows closely that of previous work on Voronoi refinement meshing p2 H2] . We 
will show that certain geometric properties are satisfied after each point is inserted by the 
algorithm. The Feature Size Invariant (Theorem pj guarantees that no two output points 
are too close to each other compared to the feature size of the input points. It is used to prove 
the termination of the algorithm and the asymptotic optimality of the output size (see |16[ Ch. 
3] for a general version of this fact). 

The other invariant we maintain is the Quality Invariant (Theorem [HI). It guarantees 
that the points in each set of the hierarchy are well-spaced. This is a critical requirement for 
a point set to be hierarchically well-spaced. It also guarantees, among other things, that the 
approximate Delaunay graph has constant bounded degree throughout the algorithm. 

These invariants will apply to the whole hierarchy of point sets. The sets in the hierarchy do 
not interact, because the approximate Delaunay graph has no edges between points in different 
sets. As a result, it suffices to prove each invariant individually, for each set in the hierarchy. 

At each step of the algorithm, there is an ordering of the inserted points. This ordering is 
mostly determined by the order in which the points were added. The only exception is that if 
an input point p causes a nearby Steiner point q to snap, then p is given the index of q in the 
updated ordering. 

Let Mi, . . . , Mm be the sequence of ordered point sets for each step of the algorithm. Simi- 
larly, let Pi = Mi n P for each i = 1 . . . N . Each step either adds an input point, adds a Steiner 
point, or snaps a Steiner point to an input point. If the ztti step is one of the first two cases, 
then Mi is formed from Mj_i by appending a new point to the end of the ordering. If there is 
a snap then Mi is formed from M,-_i by swapping one Steiner point for an input point. It is 
clear that for all i < j, Pi C Pj. 

The reader familiar with other analyses of Delaunay or Voronoi refinement schemes will 
notice a difference here in how we define the steps of the algorithm. In previous work, one 
considered a single ordered point set and reasoned about the prefixes of that ordering. This is 
also standard in analyses of other incremental geometric algorithms. In our case, we may undo 
some Steiner point insertion if it was too close to an input point, so we need a more nuanced 
analysis. The payoff is that we avoid the preemptive point location for Steiner points used in 
sparse Voronoi refinement [9], and we get by without any extra point location data structures. 

4.1 The Feature Size Invariant 

We want to show that at each step in the algorithm, the feature size with respect to the 
inserted points is larger than some constant times the feature size with respect to the input 
points. Recall that the X(v,M) is the predecessor radius of the point v in the ordered set 
M. Usually, \(v, Mi) = X(v,Mj) for any Mj and Mj containing v. However, for some points, 
snapping a Steiner point could cause the predecessor radius to change. If it goes up, that only 
makes the feature size invariant easier to satisfy. If it goes down, we need to carefully check 
that it does not go down by too much. The following standard lemma shows that bounding the 
predecessor radius suffices to bound the feature size. 



Lemma 2. Let M C M. d be an ordered set with order relation -<. Let f : M. d — > R 6e any 
1-Lipschitz function and let K > 6e any constant. If for all v G M, /(u) < (if — l)A(w, M), 
then for all v G V, f(v) < Kf M (v). 

Proof. Fix any v G M and let -u = NN(u, M). If u -< v then ffcf(u) = d(u,i>) = X(v,M), and 
so, /(«) <(K- l)X(v, M) = {K - l)hi{v) < Kf M {v). If v -< u, then X(u, M) < d(u, v), and 
so, f(v) < f{u) + d(u, v) <(K - l)A(u, M) + d(u, v) < Kd(u, v) < Kf M (v). D 

Theorem 3 (The Feature Size Invariant). For all i = 1 . . . N , and for all v G Mi, fp^v) < 
Kf M .(v), where K = ^. 

Proof. Lemma [2] implies that it will suffice to bound the predecessor radius of each point in each 
ordering. Define the subset Xi C Mi to be the points for which some previous snap caused a 
decrease in the predecessor radius, i.e. Xi := {v G Mj : X(v, Mj) < X(v, Mj-i) for some j < %}. 
We will proceed by induction. The inductive hypothesis is that for all k < i and all v G M^, 

f K P X(v,M k ) ifv£PiUX k 
P * [V > ~ \ K s X(v, M k ) iiveM k \P t \X k 

where K$ = ^| and Kp = ^r^- Note, this choice of Kp and K$ implies 

2K 
K P = K - 1, K s = + 1, and K P = 2K S + 1. (1) 

r 

We will now prove that this hypothesis also holds for Mj. There are three cases to consider, 
a new Steiner point, a new input point whose nearest neighbor is also an input point, or a 
Steiner point snapped to an input point. 

Case 1: Adding a Steiner point If we are adding a Steiner point v, then the inductive 
hypothesis continues to hold for all points in Mj other than v. We need only show that it holds 
for v as well. Let x be the vertex whose cell had bad aspect ratio. Since v is in the Voronoi cell 
of x before insertion, X(v,Mi) = d(x,v). So, we may assume 2d(x,u)/fM i _ 1 (x) > r. Moreover, 
by induction and Lemma [2l fp,(x) < KiMi^ x {x) = Ki^j^x), where the last equality follows 
because v is not the nearest neighbor of x in Mj. So, we have that fp^x) < 2Kd(x, v)/t. Using 
these facts and the Lipschitz property of fp i , we derive the following. 

f Pi (v)<f Pi (x) + d(x,v) 

< \-^r + M d ( x ' v ) = K s X(v,Mi). 

Case 2: Adding an input point (no snapping) Since the only change in the ordering 
is to append the new vertex v G Pi to the end of the ordering, the inductive hypothesis will 
continue to hold for all points other than v in Mj. Let u = NN(u, Mj). Since we did not snap, 
it must be that u G Pj. So, fp^v) < d(u, v) = X(v, Mi) < KpX(v, Mi). So, the hypothesis holds 
for v as well. 

Case 3: Snapping a Steiner point to an input point Let v be the new point added 
and let x be the Steiner point that was removed. It will suffice to check that the inductive 
hypothesis holds for v as well as any points for which the predecessor radius has been decreased 
by the snap. For such a point u, it must be that x precedes u in Mj_i and so using the triangle 
inequality and that d(x,t>) < d(n, v), 

X(u, Mj_i) < d(x, u) < 2d(n, v) = 2X(u, Mi). (2) 



10 



Case 3(a): u = v By the definition of x, we have the d(x, v) < X(v, Mi). So, by the triangle 
inequality, A(x, Mj_i) < X(v,Mi) + d(x,v) < 2\(y,Mi). If x G AQ-i, then there is some point 
y € Pi such that A(x,Mj_i) = d(x,y). Using Q, we get the following bound. 

f») < d(v,y) < d(v,x)+d(x,y) < 3A(t>,M) < K P X(v,Mi). 

Otherwise, if x g Xj_i, then by induction and Q, fp t (x) < K s X(x,Mi_i) < 2K s X(v,Mi). So, 
we get that f P .(» < f P . (x) + d(x, v) < (2K S + l)X(v,Mi) = K P X(v,Mi). 

Case 3(b): u G P t \ {v} f Pi («) < d(«,«) = A(tt, Afc) < K P X{u,Mi). 

Case 3(c): u G Xj_i \ Pi The predecessor radius is determined by an input point for all 
points of Xi, so fp 4 («) < A(ti,Mj_i). Using Q, we have fpj('u) < A(ti,Mj_i) < 2A(n,Mj) < 
K P X(u,Mi). 

Case 3(d): u G X { \ Xj_ x \ P,. By induction and (|2|, we have fp.(«) < K 5 A(n,Mi_i) < 
2 J FC s A(u, M^ < i^ P A(u, M^. □ 

There is one special case of the preceding theorem that is sufficiently useful that we state it 
below as its own lemma. Its proof is the same as above, but only requires the first case. 

Lemma 4. Let Mi, . . . , M^ be a contiguous sequence of ordered sets produced by the algorithm in 
which each was formed from the previous one by adding a Steiner point. Then, for all j = i . . .k, 
and for all v G Mj, fAf, ; (w) < K^M {v), where K = ^£4. 

4.2 The Quality Invariant 

We now show that at all times, i.e. for each M±, . . . M^v, the points in each set of the hierarchy 
are well-spaced for a quality constant that does not depend on the dimension. In particular, 
this means that each Mj is indeed a hierarchically well-spaced point set. 

Theorem 5 (The Quality Invariant). There exists a constant r" independent of ' n and d, such 
that at all times, the set of inserted points is t" -well-spaced. 

The proof of this fact follows exactly the pattern of previous work [9j [12] and is broken into 
Lemmas [6] and [8] below. First we show that inserting an input point can decrease the aspect 
ratio from r to a constant r'. Then, during the refinement the aspect ratio never dips below r" . 
By construction, the refinement phase will end with aspect ratio at most r as computed using 
the approximate farthest corner routine. Since this is approximation, the true aspect ratio of 
the cells may be worse by a small constant factor. 

Lemma 6. Adding one point to a r-well-spaced set using the Insert routine results in a r'- 
well-spaced set, where t' = 2Kt. 

Proof. Let p be the newly inserted point. The only way for the aspect ratio of a point q to 
increase is if its inradius decreases. This can only happen if the nearest neighbor is p. If p 
was added to the Voronoi cell of q, then q must be an input point for otherwise, the algorithm 
would have snapped q to p. Moreover, for q G P, we only insert p if d(p, q) > r(q)/K, so r{q) 
decreases by at most a factor of 2K. For any other q whose Voronoi cell's inradius decreases, p 
cannot be more than twice as close to q than the previous nearest neighbor of q. In that case 
the aspect ratio of q goes down by at most a factor of 2. □ 

To prove that a sequence of Steiner points added during refinement maintains a bound on 
the quality, we will use the following lemma from Hudson et al. [9l Lemma 6.1]. 
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Lemma 7. Let M be a set of r' -well-spaced points, let c be a point in the convex closure of M, 
and let B = ball(c, r) be a ball that does not contain any points of M in its interior. Then, for 
allx £ B, f A fO) > j^m- 

Lemma 8. Let Mi and M^ be the points of a mesh before and after a sequence of Steiner points 
were added, where Mi is t' -well-spaced points. Then, for each j = i . . .k, Mi is t" -well-spaced, 
where t" = 32Kt' 2 . 

Proof. Fix any intermediate mesh Mj and any vertex v G Mj. Let R = R(v) be the outradius 
of Vor Mj(v)- So, the aspect ratio of ~Votmj(v) is bounded using Lemmas El and [7| as follows. 

2R . 2KR . „„„ /2 



aspect M (v) = < < Z2Kt . 

3 ImAv) l Mi {v) 



D 



A Note About Constants This section is optimized for concise proofs rather than optimal 
constants. The main slack that can be tightened by a more detailed analysis comes from the 
constant r' defined in Lemma [6| Although the constant is tight, at most two points can achieve 
the bound at any given time. In the application of Lemma [7j it is assumed that all points may 
have aspect ratio r', and this assumption leads to a much looser bound than could be realized. 

5 Correctness and Running Time 

We now show that the described algorithm correctly computes a hierarchically well-spaced point 
set and maintains a superset of the edges in the Delaunay graph in each set in the hierarchy. 

5.1 Delaunay Neighbors and Degrees 

In this section, we prove the two most important properties of the approximate Delaunay graph, 
namely that the graph contains the true Delaunay graph and also has bounded degree. We will 
assume that the points are in sufficiently general position that the Delaunay triangulation is 
indeed a triangulation. If this is not the case, any infinitesimal (or symbolic) perturbation will 
make it so. 

Theorem 9. The approximate Delaunay graph produced by our algorithm contains all edges of 
the Delaunay triangulation of each point set in the output hierarchy of well-spaced points. 

Proof. The proof is by induction on the number of insertions and deletions performed by the 
algorithm. We show that at each step, the approximate Delaunay graph contains all of the 
Delaunay edges of each set in the hierarchy. 

As a base case, note that the starting cage connecting all pairs of vertices clearly contains 
all Delaunay edges. Assume inductively that the set Mj_i from the first i — 1 insertions and 
deletions satisfies the inductive hypothesis. 

Let Mi be the set formed by inserting p into M,_i. By the Quality Invariant (Theorem pi), 
there are no Delaunay edges (p, q) longer than mm{2 R(p), 2 R(q),2T"r(p)}. So, when the algo- 
rithm prunes such edges from the approximate Delaunay graph, it does not remove any Delaunay 
edges. It will therefore suffice to prove that all new Delaunay edges formed by inserting p are 
added to the approximate Delaunay graph. All such edges are necessarily incident to p. At 
least one such edge is the edge from p to its nearest neighbor q = NN(p, Mj_i). Note that the 
neighbors of p in Del^ induce a connected subgraph of the Delaunay triangulation. Indeed, 
by Delaunay/Voronoi duality, this graph is isomorphic to the edges of the polytope polar to 
VoTMi(p)i arid polytopal graphs are connected. So, the graph search starting at q will find all 
of the Delaunay neighbors of p. 
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Now, we consider the case of a point deletion that occurs when we snap. Let Mj be the 
set of points formed from removing a point p from Mi—\. In this case, the only new Delaunay 
edges in DelM; that are not in ~DelM i _ 1 are those that connect two neighbors of p. Since we 
tentatively add all such pairs before pruning those that cannot be Delaunay, the updated graph 
again contains all Delaunay edges. □ 

We now prove an upper bound on the degree of a node in the approximate Delaunay graph 
at any time in the algorithm. 

Lemma 10. In the approximate Delaunay graph of a well-spaced mesh with quality t" , each 
point has at most (3r" 2 ) neighbors. 

Proof. Let p any vertex in the approximate Delaunay graph. Let t be its degree and let q\ , . . . , q£ 
be its neighbors. Consider the inballs of q%, . . . , q#. As an immediate consequence of the pruning 
condition, each neighbor qi has inradius r(%) < |d(p, qi) < T"r(p). Similarly, r(qi) > r(p)/r". 
Moreover, 

d(p,qi) < 2r"r(p). 

Therefore, all the inballs lie within a ball of radius d(p, qi) +maxj{r(gj)} < 3r"r(p) centered at 
p. Since the inballs are disjoint, the following volume packing argument establishes the claim. 

vol(ball(3T"r(p))) = „ 2 d 
~ vol (ball(r(p)/r")) l ' ' 

□ 

5.2 Calculating Approximate Aspect Ratios 

We now prove the correctness of our method of approximating the aspect ratio of a Voronoi cell 
by finding an approximate farthest corner. 

First, we show that we can approximate the farthest corner of a Voronoi cell in any specified 
direction in the desired running time by using the ellipsoid method. For each neighbor q of p, 
we create a linear constraint that x £ M lies on the same side of T-L{p,q) as p (here, Ti(x,y) 



denotes the hyperplane bisecting xy). Let C be the set of all such constraints. By Lemma 10 
there are 2 °( > constraints. Moreover, it is clear that the feasible region of the linear program 
defined by these constraints is precisely Vor(p). 

Theorem 11. Let p be a point in the mesh which is, without loss of generality, located at the 
origin, and let c G M d be a vector with \\c\\ = 1. Suppose w is the farthest corner of\oi{p) in the 

direction of c, i.e., w is the point ofVor(p) that maximizes w T c. Then, in 2 ( 'O (log ( — ) ) 

time, we can calculate a point z £ Vor(p) such that z T c > (1 — e)w T c. 

Before we prove the above theorem, we prove a lemma. 

Lemma 12. Let p, c, w, and C be as before, and assume s = w c. Write r = r(p) and 
R = R{p). Suppose that t £ [r, rr"] and s > j^r. Then, in 2 °( 'O (log ( ^r ) ) time, we can 
produce output according to the following rules: 

1. If t < s (l — |), then the output should be a point x £ R that satisfies the constraint set 
C U {x T c > t}. 

2. Ift> s, then the output should indicate "NOT FOUND." 

3. If s (l — |j < t < s, then the output can be either a point x £ R satisfying C U {x c > t} 
or "NOT FOUND." 
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(a) The shaded volume shows the cone 
formed by joining w to ball(p, r). It is 
contained in the Voronoi cell of p. 




(b) The shaded volume shows the cone 
formed by taking the intersection of the 
cone in (a) with the half space {x £ M. d : 
x T c > t}. The dotted arrow indicates the 
direction given by c. 



Figure 1 : Bounding the volume of the feasible region during the ellipsoid method 



Proof. The idea is to run the ellipsoid method for LP feasibility, but we make sure to terminate 
the process early after 0(d 2 log(r"/e)) iterations. At this point, if the center of the final ellipsoid 
is feasible, we output the center; otherwise, we output "NOT FOUND." It is evident that rules 
(2) and (3) will be satisfied by the algorithm. 

Let us show that (1) is also satisfied. Suppose t < s (l — |). Then, the corresponding 
LP is feasible. We shall now bound the volume of the feasible region TL from below. Recall 
that the ball ball(p, r) of radius r centered at p is entirely contained in Vor(p). Thus, the d- 



dimensional cone joining apex w to ball(p, r) is contained in Vor(p) (see Figure la). Hence, the 
intersection of this cone with {i £ R : x c > t} (which is itself a cone with an ellipsoid base) 



is contained inside TZ (see Figure lb). Using elementary geometry, we find that the volume of 

this intersection is at least 

1„ r d - l {s-t) d 



d ( s 2 



2^— ' 

p 2 



where Cj is the volume of a j-dimensional unit ball (see Appendix |B| for a detailed derivation) 
Using the assumptions t < s (l 
and is at least 



|) and s > T^r, we see that the above quantity is defined 



d Cd ~ l \2 



„d— l„d 



'-)' 



The above quantity is minimized for s 
region is bounded from below by 



y/dr. Thus, it follows that the volume of the feasible 




Moreover, by Theorem p3l we know that ball(p, rr") contains Vor(p) and, therefore, the feasible 
region as well. Hence, by a standard analysis of the ellipsoid method (see [7]), the number of 

iterations required to correctly output a feasible point is O ( d log f voiiji)' ) ) > wmcn we 
compute as 



O dlog 



C d (rr 



ll\d 



l°d- 



(er/2) c 



O d 2 log 
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Terminating the ellipsoid method after O I d 2 log ( %- ) J iterations therefore suffices to obtain 
the desired output, hence establishing condition (1). 

Finally, since performing each iteration of the ellipsoid method can be done in 2°W time, 
the total running time is 2°^0 (log ( — ) ) , as desired. □ 



Now, we are ready for the proof of Theorem 11 



Proof of Theorem 11. We show how to obtain z. Let s = w c. The trick is to run the feasibility 
algorithm of Lemma 12 on the linear program with constraint set C U {x T c > t} for different 



values of t. Let L 



2(r"-l) 



Then, let to, ti, t2, ■ ■ ■ , ti, G [r, rr"] be defined by U = r + 



L 



In particular, we do a binary search to find i £ {0, 1, . . . , L] such that the algorithm of Lemma 12 
returns a solution for t = U but reports "NOT FOUND" for t = U + \. Then, we simply output 
the z corresponding to the obtained solution for t = U. If such an i does not exist (i.e., the 
algorithm does not return a feasible solution for any i), we simply output z = re. Note that 

the total runtime of this overall procedure is 2 ' >0 ( log ( -^r ) ) • 

It remains to show that z c > (1 — e)s. 

1. If s < r^r, then note that any z that is outputted satisfies z c > r > (1 — e)s. 

2. Suppose s > i^r. Let i be the index returned by the binary search for which t = U 
returns a solution but t = tj + i returns "NOT FOUND." Clearly, we must have that 
U+l > (l — |) s. Thus, we get the desired bound on ti as follows. 

rr" — r ( e \ er 

t,=t m -^— >(l--j,-- 

> U-o )r-- = (l-e)r. 



D 

Now, the farthest corner in the direction of an arbitrary unit vector c is, in general, not the 
farthest corner of the entire Voronoi cell. Hence, we need the following lemma. 

Lemma 13. Let p be a point in the mesh with q as the farthest corner of its Voronoi cell. Let 
< 6 < 7r/2 and let c G S such that the angle between c and q — p is < 6. Then, applying 



the underlying algorithm of Theorem 11 with this choice of c produces a point z inside Vor(p) 



for which d(p, z) > (1 — e)(cos#)d(p, q). 

Proof. As before, assume without loss of generality that p is located at the origin. Let w be 



the farthest corner of Vor(p) in the direction of c. The algorithm of Theorem 11 produces a 
point z inside Vor(p) such that z T c > (1 — e)w T c. The desired claim follows since w T c > q T c > 
d(p, q)cos9. □ 

As a corollary, we obtain the following: 



Corollary 1. If q is the farthest corner ofYoi(p), then repeating the algorithm of Theorem 11 
for each c S rj— CAGE and taking the point z that is farthest from p among all outputted points 
gives us d(p, z) > (1 — e) ( 1 — \ ) d(p, q). 



Proof. Consider the point z 1 = (z — p)/d(p, z). Since z 1 € S d 1 , there exists c E ry— CAGE such 
that d(z',c) < r). Since z' and c are both unit vectors, elementary geometry imp 



ies that the 
implies the 
desired claim. □ 



angle 9 between z' and c satisfies cos 6 = 1 z ^ ^ 1 — \- Thus, Lemma 
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5.3 The Greedy Walk has Constant Length 

The greedy walk is named as such because it uses the predecessor pairing from the greedy 
permutation to start the walk. This is in contrast to previous walking-based point location 
schemes that use randomization. Also, it walks through the Voronoi diagram rather than the 
Delaunay simplices and therefore needs no hyperplane tests. In this section, we show that the 
greedy walk only takes a constant number of steps before terminating. 

Lemma 14. For an input point p, the greedy walk from 4>{p) to p only takes 2 °( > steps. 

Proof. Let M denote the current set of inserted points at the time we insert p. Let S be 
the line segment p<j)(p), and let £ denote its length. Let m be the midpoint of S. The proof 
will be by a volume packing argument. It will suffice to show that any Voronoi cell Vor m{q) 
intersecting S has an inradius in the range [y£, £] for some constant 7 = 1/(4^(1 + r)) and also 
that d(q, m) < 2£. Thus, the inballs are contained in ball(m, 3^) and their number is at most 
(3/ 7 ) d = 2°( d ). 

First we observe that for any x £ S the feature size fp t (x) is at least 1/2. Otherwise, there 
must be a point y £ Pi-i that is closer than £ to either pi or (j)(pi), but this would violate the 
greedy ordering. 

Let q be any point of M such that Votm(q) intersects S. So, some x £ S is closer to q 
than to 4>{pi) and by the triangle inequality, d(q, 4>(pi)) < 2£. The inradius r{q) of VorM(^) is 
at most half the distance from q to (j)(pi) and thus is at most £. To lower bound r(q) we use 
Theorem [3] and the quality of M as follows. 

Tin) = \h4(q) > ^pM > ^ f PM - d (^) 
1 

Thus we get that r(q) > £/{AK{1+t)) = j£. The result follows from the aforementioned volume 
packing argument. □ 

5.4 Running Time 

Theorem 15. The running time of our algorithm is 2°( d '(nlogn + m). 

Proof. The expected running time of the greedy permutation algorithm is 2°(^nlogn. Each 
step of the algorithm adds a new point so it suffices to show that each insertion takes only 



20(d) tim e _ Lemmas 14 and 10 imply that the point location work is only 2°^ d > per input 



point. Moreover, the degree bound in Lemma 10 also implies that there are only 2 ' > vertices 



to update with any insertion. Each takes only 2°^ d > time to locate the new neighbors and 



recompute the approximate farthest corners, by Theorem 11 Thus, the total running time 



after preprocessing is 2°( d 'm. □ 

6 Conclusion 

We have presented several new ideas for mesh generation. The most relevant are the use of 
linear programming to find large empty regions in the mesh without resorting to enumerating 
the simplices of the Delaunay triangulation and the use of greedy permutations as a preprocess 
for walking-based point location. The latter idea alone could replace the rather complicated 
range net methods of previous work to get running times independent of the spread [12J. 
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A Pseudocode 

We present the pseudocode for the algorithm with the goal of giving the most concise description 
that both allows a rigorous analysis and could guide an implementation. The blocks of code are 
whitespace delimited, i.e. indentation indicates scope. 

We assume a graph data structure supporting dynamic insertions and deletions using the 
methods InsertEdge, DeleteEdge, Insert Vertex, and DeleteVertex. The vertices are 
the geometric points. The approximate Delaunay neighbors of a vertex v can be accessed from 
Nbr(u). 

The main routine for constructing a r-well-spaced set of points, given an input set P and 
quality parameter r, is as follows. We indicate r as an input parameter, but we assume that it 
is available to all of the methods without explicitly passing it. 

WellSpacedPoints (P: points, r: quality parameter) 
(P, 4>) <— ApproximateGreedyPermutation(P) 
Insert Vertex (pi ) 
NewLayer(p 2 ,pi) 
for i = 3 to \P\ 

Insert (pj) 

RefineQ 

The following is a routine for considers a new point p for insertion. 

Insert (p: point) 

q «— Greedy Walk (p, 4>(p)) 

if q is a Steiner point then Snap(p, q) 

else if d(p,q) < hr{q) then NewLayer(p, q) 

else RegularInsert(p, q) 

RegularInsert inserts a new point p that lies in the Voronoi cell of an already-inserted point 
q. Snap snaps a Steiner point q to an input point p which is then inserted. NewLayer adds 
a new layer to the hierarchy, places point q inside, and then inserts the input point p. 

RegularInsert (p: point, q: point) 
Insert Vertex(p) 
r{p) *-d(p,q)/2 
Do a graph search starting from q: 

At a vertex v, if d(v,p) < 2mm{R(v),T"r(p)} 

then InsertEdge (p, v) and search the neighbors of v 
UpdateAspect(p) 
PruneEdges(p) 
for each v E Nbr(jj) 

UPDATEASPECT(t;) 

PruneEdges(u) 
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Snap (j>: point, q: point) 

q' <- argmin p , eNBR(g) d(p,p') 
DELETE(g) 

RegularInsert(p, q') 

NewLayer (p: point, q: point) 

Create a new copy of 77— CAGE with center q and radius 2d(p, q) 
for each pair of vertices (u, v) 

InsertEdge(u, v) 
RegularInsert(p, q) 

Next is a procedure to remove a point p that has previously been inserted. The routine is used 
during snapping, when a Steiner point must be removed to make way for a nearby input point. 

Delete (p: point) 

N <r- NBR(p) 

Delete Vertex (p) 
for each u,v G N 

InsertEdge(u, v) 
for each v G N 

UpdateAspect (-u) 
for each t)6JV 

PRUNEEDGES(f) 

The Refine procedure repeatedly adds Steiner points at the approximate farthest corners of 
any cells whose aspect ratio is too large. 

Refine () 

while Q re fi ne is nonempty 

q <- pOp(Q r efinc) 

if R(q)/r(q) > r then RegularInsert(ApproximateFarthestCorner((/), q) 

The next two methods update the data structure. UpdateAspect computes the inradius and 
(approximate) outradius of a vertex, and it puts that vertex on the refinement queue if its aspect 
ratio is greater than r. 

UpdateAspect (p: point) 
a ■<— R(p)/r(p) 
r(p) i- \ min ggNBR(p) d(p, q) 
x <(— ApproximateFarthestCorner(p) 
R(p) <— d(x,p) 
if R{p)/r{p) > r > a then push p to Q rc fi nc 

PruneEdges (p: point) 
for each q G Nbr(j>) 

if d(p, q) > 2 min{R(p) , R(q)} then DeleteEdge(p, q) 

The next procedure does a greedy walk along p'p, starting at p', and returns the point whose 
Voronoi cell contains p. It should be noted that the notation T~L(u,v) denotes the hyperplane 
that bisects uv. 
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Greedy Walk (p: point, p': point) 

^current > P 
V next «- p' 
W 4— p* 

while ttj lies on p/p 
s 4— oo 
for each q' / f currcil t in NBR(u next ) 

if Wp intersects 7i(v ncx t, q') and (p — w,q' — v ncxt ) > then 
t' 4-wpCM-L{v ncxt ,q') 
if d(w, t') < s then 
t 4- t' 
qi-q' 
s <- d(w,t') 
if s ^ oo 
u> <— £ 

^current ^ ^next 
« next <- Q 

else 

return f currcnt 

The following routine takes in a point p and locates an approximate farthest corner of Vor(p), 
i.e., returns a point z inside Vor(P) such that d(p, z) is within a factor (1 — e) f 1 — \ J of the 
distance from p to the true farthest corner: 

ApproximateFarthestCorner (p: point) 
r ■<— r(p) 
z ■<— p 
for each c G r?— CAGE 



w -^ re 



L 



2(t"-1) 



fc^L + 1 
while k — j > 1 



TT 



t ^— r + , 

Consider the LP feasibility problem: 

xeR" 1 subject to: V<? G NBR(p), (q - p) T x < (9 ~ p) ^ (g+p) , 
c T (x — p) >t 
Run 0((i 2 log(r // /e)) iterations of the ellipsoid method 
if a feasible solution x is produced then 

J «- « 

else 

/c <(— i 
if d(p, w) > d(p, z) then z <— w 
return z 

Finally, ApproximateGreedyPermutation is a blackbox that takes in a set P of points and 
returns a permutation p±,P2, ■ ■ ■ iP\P\ of the points along with the predecessor pairing 0. 
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B Ellipsoid Volume Calculation 



In the proof of Lemma 12 we state (without proof) that the volume of the intersection of the 
cone in Figure lb is at least 

1 r d - l (s-t) d 

^1 d ~ 1 4=1' 

We shall prove that claim here. 

Without loss of generality, suppose the apex w of the cone is located at the origin and that 
the coordinates of p are (a, 0, 0, . . . , 0, s). Then, we wish to calculate the volume of the region 

d-2 

given by the intersection of the cone with {{x\,X2, • • • , xj) G M d : < x d < s — t}. Note that 
this region itself is a cone with an ellipsoid base and height s — t. 

Let us determine the ellipsoid base. Note that the boundary of the ellipsoid is given by the 
intersection of the boundary of the cone with the hyperplane x d = s — t. Note that the angle 

formed by the side of the cone with the axis is 9, where cos 9 = \/ % s ,~i r ■ Hence, the equation 
of the cone is 



(xi,x 2 ,..-,x d ) ) = \Jx\ 



a - " - S l x '2 i ^.2 i i ^,2 



or. 



,— ^,0,0,. ..,0, . ,{xi,x 2 ,..-,x d ) ) = Jxf + x^ \-xjcos9, 

V« + s s v ' vu + s z 

d-2 



ax l sx d \ / 2 , 2 , ,2\ 2/i ( 2 i . 2\ a + s ' r 

+ ; = (x 1 +x 2 + --- + x d )cos 9 = {x l + --- + x d )- 



V^T^ v^+TV Vi * a ' Vi a ' a 2 + s 2 ' 

where ( , ) denotes the standard inner product on M. d . To compute the intersection with the 
hyperplane x d = s — t, we plug in x d = s — t into the above equation to obtain 

1 s 2 —r 2 J 1 . n 9 n . 

+ ; : : : 2 (4 + 4 + • • • + x 2 d _ x ) = 1. 



r(s-t)Va 2 +s 2 -r 2 \ 2 ( r(s-t) x ' 



Hence, the ellipsoid base of the cone whose volume we desire has semiaxes with lengths 

r(s — t)Va 2 + s 2 — r 2 r(s — t) r(s — t) r(s — t) 

s 2 -r 2 ' \/s 2 - r 2 ' \/s 2 - r 2 ' ' Vs 2 — r 2 ' 

V v ' 

d-2 

Since the height of the cone is s — t, it follows that the volume is 

1. .„ (r(s - t)Va 2 + s 2 - r 2 \ ( r{s - t) \ d ~ 2 1 r dr - 1 (s-t) d 

d {s ~ t)c " { ^7 2 ) KW^T 2 ) ~ ^- 1 (s2 _ r2) ^- 
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